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We report on localized patches of cellular hexagons observed on the surface of a magnetic 
fluid in a vertical magnetic field. These patches are spontaneously generated by jumping 
into the neighborhood of the unstable branch of the domain covering hexagons of the 
Rosensweig instability upon which the patches equilibrate and stabilise. They are found to 
co-exist in intervals of the applied magnetic field strength parameter around this branch. 
We formulate a general energy functional for the system and a corresponding Hamiltonian 
that provides a pattern selection principle allowing us to compute Maxwell points (where 
the energy of a single hexagon cell lies in the same Hamiltonian level set as the flat 
state) for general magnetic permeabilities. Using numerical continuation techniques we 
investigate the existence of localized hexagons in the Young-Laplace equation coupled 
to the Maxwell equations. We find cellular hexagons possess a Maxwell point providing 
an energetic explanation for the multitude of measured hexagon patches. Furthermore, 
it is found that planar hexagon fronts and hexagon patches undergo homoclinic snaking 
corroborating the experimentally detected intervals. Besides making a contribution to 
the specific area of ferrofluids, our work paves the ground for a deeper understanding of 
homoclinic snaking of 2D localized patches of cellular patterns in many physical systems. 


1. Introduction 

Spatial localization is proving to be key to understanding various coherent features 
seen in fluids comprising Taylor-Couette flow (Abshagen et al. 2010), plane Couette 
flow (Schneider et al. 2010; Chantry & Kerswell 2014), binary fluid convection (Moses 
et al. 1987; Batiste et al. 2006; Mercader et al. 2013), and transition-to-turbulence (Pringle 
et al. 2014). See also Knobloch (2008, 2015), Dawes (2010) and Descalzi et al. (2011) for 
broader reviews. In the Homoclinic Snaking scenario (Pomeau 1986; Woods & Champ- 
neys 1999; Coullet et al. 2000) each alternating turn of a “snake” in control parameter- 
phase space is correlated with the emergence of a further interior cell of the localized 
pattern. This scenario has proven to be an important mechanism for localization in more 
than hundred theoretical studies. However, experimental evidence is limited, so far, to 
buckling (Thompson 2015), a semiconductor laser (Barbay et al. 2008), vertically vibrated 
media (Umbanhowar et al. 1996), nonlinear optics (Akhmediev & Ankiewicz 2005), gas 
discharges (Purwins et al. 2010), a light valve (Haudin et al. 2011) and in the context of 
fluids, the Taylor-Couette flow (Schneider et al. 2010; Abshagen et al. 2010). A pending 
problem is a comparison of experiment and theory of homoclinic snaking for fully local¬ 
ized 2D patterns (not just structures that are localized in a single spatial direction) in a 
physical fluid system. 
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Figure 1. A patch of seven ferrofluidic spikes seen together with the upper coil (a), and zoom 
in high density range (b). Photographs courtesy of Robin Maretzki (a) and Achim Beetz (b). 


A step towards filling this gap is investigating the surface instability of an electrically 
or magnetically polarizable fluid, because spatial inhomogeneities and noise are very low. 
When a static layer of magnetic fluid is subjected to a vertical magnetic field, it becomes 
unstable when a critical value Be of the magnetic induction is surpassed. Due to this 
Rosensweig instability (Cowley & Rosensweig 1967) a regular hexagonal pattern of cellu¬ 
lar spikes emerges. By applying local magnetic pulses, localized radially symmetric spikes 
could be generated on the free-surface (Richter & Barashenkov 2005). In this paper we 
demonstrate that various localized hexagon patches like in figure 1 emerge spontaneously 
if the homogeneous field is switched into a regime where we expect to find homoclinic 
snaking. So far the modeling of localized spikes is limited to finite element simulations 
(Lavrova et al. 2008; Cao & Ding 2014) and a simple reaction-diffusion model (Richter & 
Barashenkov 2005). However, the Rosensweig instability involves a free-surface and so the 
justification of reaction-diffusion theory is problematic. Solving the full dynamic equa¬ 
tions, we present numerical evidence that the localized hexagon patches undergo homo¬ 
clinic snaking similar to that observed in reaction-diffusion theory (Woods & Champneys 
1999; Coullet et al 2000). 

The article is outlined as follows. In §2, we describe the experimental setup, which is 
followed by a characterization of the ferrofluid (§3) and the experimental results (§4). 
We then introduce in §5 the mathematical model for the experiment that we numerically 
solve and present in §6 the numerical results. Section 7 then compares experimental and 
theoretical results. We finish with conclusions and outlook (§8). 


2. Experimental setup 

Our experimental setup is sketched in figure 2(a). An X-ray point source emits radiation 
vertically from above through the container with the ferrofluid, which is placed midway 
between a Helmholtz pair of coils. Underneath the container, an X-ray camera records 
the radiation passing through the layer of ferrofluid. The intensity at each pixel of the 
detector is directly related to the height of the fluid above that pixel. Therefore, the 
full surface topography can be reconstructed after calibration (Richter & Biasing 2001; 
Gollwitzer et al 2007). 

The container, which holds the ferrofluid sample, is depicted in figure 2(6). It is a 
regular octagon machined from aluminium with a side length of 77 mm and two concentric 
inner bores with a diameter of 140 mm. These circular holes are carved from above and 
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Figure 2. Setup of the apparatus for measuring the Rosensweig instability. Scheme of the 
assembled setup (a), and exploded view of the container with the coils (b). 


below, leaving only a thin base in the middle of the vessel with a thickness of 2 mm. 
On top of the octagon, a circular lid of aluminium is placed, which closes the hole from 
above (see figure 2b). Each side of the octagon is equipped with a thermoelectric element 
QC-127-1.4-8.5MS from Quick-Ohm. These are powered by a 1.2 kW Kepco KLP-20-120 
power supply. The hot side of the thermoelectric elements is connected to water cooled 
heat exchangers. The temperature is measured at the bottom of the aluminium container 
with a PtlOO resistor. The temperature difference between the center and the edge of 
the bottom plate does not exceed 0.1 K at temperature 10.0 °C measured at the edge of 
the vessel. A closed loop control, realized using a computer and programmable interface 
devices, holds the temperature 0 of the vessel constant with an accuracy of 10 mK. 

The container is surrounded by a Helmholtz pair of coils, thermally isolated from the 
vessel with a ring made from a flame resistant composite material (FR-2). As shown 
in figure 2 the diameter of the coils is intentionally smaller than the diameter of the 
vessel in order to introduce a magnetic ramp. With this arrangement the magnetic field 
strength falls off towards the border of the vessel, where it reaches 80% of its value 
in the center. The ramp is introduced in order to minimize the differences between the 
amplitudes obtained in a finite and an infinite geometry. For details of the magnetic 
ramp see Knieling et al. (2010). Filling the container to a height of 5 mm with ferrofluid 
enhances the magnetic induction in comparison with the empty coils for the same current 
/. Therefore B{I) is measured immediately beneath the bottom of the container, at the 
central position, and serves as the control parameter in the following. 


3. Properties of the Magnetic Fluid 

For the experiments we use the ferrofluid APGE32 from Ferrotec Co. Table 1 lists 
the important material parameters. The density p was measured utilizing an electronic 
density meter (DMA 4100) from Anton Paar. The surface tension was measured by 
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Quantity 


Value 

Error 

Unit 

Density 

P 

1168.0 

±1 

kgm“^ 

Surface tension 

a 

30.9 

±5 

mNm“^ 

Viscosity at 10 °C 

g 

4.48 

±0.1 

Pas 

Initial permeability from (3.2) 

Pr,i 

4.57 

±0.005 


Saturation magnetization from (3.2) 

Ms 

23.15 

±0.8 

kAm-^ 

Critical induction from (3.6) 

Rc,oo 

10.83 

±0.1 

mT 

Critical induction from experiment 

Rc,exp 

11.23 

±0.1 

mT 

Table 1. Properties of the magnetic fluid APGE32 (Lot G090707A) 

from Ferrotec Co. 


means of a ring tensiometer (Lauda TE 1) and corroborated by the pendant drop method 
(Dataphysics OCA 20). 

The viscosity g has been measured in a temperature range of—5 °C ^ 0 ^ 20 °C, 
using a commercial rheometer (MCR301, Anton Paar) with a shear cell featuring a cone- 
plate geometry. At room temperature, the magnetic fluid with a viscosity of = 2 Pa • s 
is 2000 times more viscous than water. The value of g can be increased by a factor of 9 
when the liquid is cooled to —5 °C. The temperature dependent viscosity data can be 
well described by the Vogel-Fulcher law (Rault 2000) 


g = go exp 


O-Oo 


(3.1) 


with go = O.48mPa'S,'0 = 1074K, and Oo = —107.5 °C, as reported in detail by Goll- 
witzer et al. (2010). For the present measurements, we chose a temperature of 6 > = 10 °C, 
where the viscosity amounts to 77 = 4.48 Pa - s according to equation (3.1). 

The magnetization, M(Hi?), where Hi? is the magnetic field in the fluid, has been 
measured by means of a fiuxmetric magnetometer (Lakeshore Model 480). A spherical 
cavity with a volume of 1 ml was constructed from perspex® and filled up to the brim 
with the magnetic fluid. By means of x-rays we checked for a bubble free filling. Figure 
3(a) presents the measured data together with a fit (solid line) by equation 




(3.2) 


utilized by Frohlich (1881) and Vislovich (1990), where M = |M|, Ms is the saturation 
magnetisation, iLp is the magnitude of the applied magnetic field in the fluid and /ir,i is 
the initial magnetic permeability. 

Fitting the data with (3.2) we obtain for the initial magnetic permeability = 
4.57, and for the saturation magnetisation Ms = 23.15 kA/m. Equation (3.2) is used to 
calculate the chord permeability 


/ich = M/ Mp + 1 


(3.3) 


and the tangent permeability 

Mta = dM/dHp + 1 = --MM-^ + 1 ^ ( 3 , 4 ) 

(l + (A^r,i-l)i|) 

which are displayed in figure 3(b) by a dashed red (blue) line, respectively. The effective 
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Figure 3. (a) Measured magnetisation curve versus the field in the fluid, and (b) related chord- 
(red dashed), tangent- (blue dashed) and effective (black solid line) permeabilities. The dotted 
line marks /ieff at He- 


permeability 

/^eff = V/^ta(^F) • Mch(^F) (3.5) 

is the geometric mean of both permeabilities (Cowley & Rosensweig 1967) and marked 
in figure 3(b) by the black solid line. 

By means of /ieff(M), the material parameters a, p and the gravitational acceleration 
g = 9.81m/s^ the critical magnetization for a semi-infinite layer of ferrofinid is approxi¬ 
mated by the implicit equation (Rosensweig 1985) 

By numerically solving the implicit equation we obtain = 6.2303 kA/m, He = 
2.3872 kA/m, and Be = po{He-\- Me) = 10.83 mT. Note that this value is just an approx¬ 
imation for a semi-infinite layer of ferrofinid. Revisiting figure 3(b), we see that /ieff drops 
considerably for increasing field. At the critical field we find /ieff (77c) = 3.2 (cf. the dot¬ 
ted line). The critical induction found by fitting amplitude equations to the experimental 
data (Gollwitzer et al 2010) is about 4% higher (cf. Table 1). 

We stress that the mapping between a linear and a nonlinear M{H) by (3.5) is only 
valid in the linear regime of pattern formation, i.e. for the onset of the instability, but 
not for the fully developed nonlinear patterns. 


4. Experimental results 

In the investigated regime the plain surface and domain covering regular up-hexagonal 
(i.e. hexagons whose maximum amplitude is postive) patterns are bistable (10.95 mT ^ 
B ^ 11.23 mT), as shown by the bifurcation diagram in figure 4 (a), which stems from a 
fit to about 170 000 measured data points (Gollwitzer et al. 2010). Panel (b) displays an 
experimental realization of the domain covering hexagonal pattern observed at the upper 
branch, the infinite extension of which has the symmetry of the dihedral group Be. 

Due to the high viscosity of the magnetic fluid the formation of Rosensweig patterns 
takes 60 seconds. This allows us to vary the magnetic induction B (the control parameter) 
and the pattern amplitude A almost independently. In this way we can enter the region 
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Magnetic induction 5(niT) 


Figure 4. (a) The bifurcation diagram (black line) of this pattern is obtained by fitting equation 
(11) in Gollwitzer et al (2010) to about 170 000 measured data. Via the path indicated by the 
arrows one obtains localized patches, as shown in (c). B is kept fixed during (1) for a time delay 
of Ti, then B is switched to (2) for a time delay of T 2 before finally being switched to (3) until 
equilibrium is reached at point (3b). Panel (b)presents a reconstruction of the height field of the 
domain covering hexagons at point (la). Panel (c) gives a reconstruction of the height field of a 
three-spike patch at point (3b) in the bifurcation diagram. 


around the unstable branch of the instability, where homoclinic snaking is suspected, 
from the side, as sketched in figure 4. 

Starting from a fiat layer, B is always switched to an over critical value Bi = 11.45 mT, 
cf. position (la) in figure 4. With time the amplitude of the domain covering hexagon 
pattern increases along path (1). After ri = 60.0s this pattern is stabilized (lb). Then B 
is switched to a value B 2 = 10.74 mT before the fold of the hexagons and the magnetic 
fluid is allowed to relax from (2a) towards the flat state. Before the magnetic fluid flattens 
completely, B is switched at (2b) to a value B^ in the neighborhood (3a) of the unstable 
branch of the imperfect transcritical bifurcation (cf. equations (1,2) in Gollwitzer et al. 
2010), marked by a dashed line and shaded region in figure 4. Now the pattern relaxes 
on path (3a) to (3b) where the height of the patch finally settles down (equilibrates) to 
a value approximately 1.5 mm higher (depending on the size of the patch i.e., for larger 
patches the height tends to that for the domain-covering hexagons) than the domain cov¬ 
ering hexagons. A movie (video 2014) shows the birth of a three spike patch. Each patch 
was at least stable for 120 seconds, as recorded by the measuring program. Moreover, we 
have checked the long term stability of exemplary patches after generation and observed 
no decay for more than two hours. 

Starting at the edge of the bistable region defined by the fold of the domain-covering 
cellular hexagon spikes and increasing Bs, we obtain a sequence of localized states ranging 
from solitary spikes to patches containing multiple spikes/hexagon cells on an otherwise 
fiat fluid surface with two to 24 cells, as shown in figure 5 for r 2 = 10.0 s. The height of 
the spikes/cells of a patch are not uniform; we observe that spikes near the centre of the 
patch are higher than those on the edge of the patch. The distance between spikes within 
a patch is approximately the same as the spatial period of the domain-covering hexagon 
pattern. Remarkably, the spikes do not spread out in space as time increases despite their 
magnetic dipole-dipole repulsion, originating from the fact that the spikes are magnetized 
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Figure 5. Localized states for T 2 = 10.0s at (a) B=11.0725mT, (b) B=11.0857mT, 
(c) B=11.1120mT, (d) B=11.1648mT, (e) B=11.1780mT, (f) B=11.1911mT, 

(g) B=11.2043mT, (h) B=11.2175 mT, and (i) B=11.2307mT. 


in parallel. In this way the patches are not an ensemble of weakly interacting solitary 
spikes in the sense of Richter & Barashenkov (2005). 

Moreover it is interesting to note the mechanism by which hexagon cells are added 
to large patches as is increased. Going from panel 5(f) to (i), we see that spikes are 
added in the middle of the lower front in order to form a completed patch with dihedral 
symmetry B 2 . It is also possible to see completed patches with De-symmetry, as shown 
in Fig. 1. 

Similar patches have been uncovered in the 2D Swift-Hohenberg-Equation as a sig¬ 
nature of homoclinic snaking (cf. Lloyd et al. 2008; Escaff & Descalzi 2009; Kozyreff & 
Chapman 2013). Also the mechanism governing the growth of a patch, as observed above, 
resemble the ones described by these authors. 

As shown in figure 6, the number of spikes in a patch increases monotonically with 
^ 3 , with clearly defined plateaus. The data related with figure 5 are marked by v 
have been recorded with a delay of r 2 = 10.0 s. Reducing r 2 to 9.5 s (□), 8.5 s(0): 

7.5s (A) shifts the curves to lower magnetic inductions: Because the starting amplitudes 
are higher, a lower is sufficient to generate a localised patch. Most importantly the 

prominent plateaus at 3, 10, 12, 14,... cells exist for all these values of r 2 . Eor the very 
short delay of r 2 = 1.0 s, the curve (o) hardly exhibits plateaus (apart from the plateau 
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Figure 6. The number of spikes vs for different delays T 2 = 1.0 s (o), 7.5s (A), 8.5 s(0)? 
9.5 s (□), and 10.0 s (v) exhibits plateaus. Plateaus at 3, 10, 12, 14 and 21 spikes, where patches 
are forming complete patches, are marked by horizontal lines. 


at 27 spikes, which is the maximal size of a hexagonal pattern in the setup). That is 
because for T 2 = 1.0 s, the regular hexagonal pattern (see (lb) in figure 4) has hardly 
decayed, and we can not sensitively explore the phase space where homoclinic snaking 
occurs. 


5. Governing equations 

In order to elucidate the snaking process we numerically solve the full 3D dynamical 
equations. The field in the fluid is modeled by Maxwell’s equations with no free currents 
and a Young-Laplace equation for the free-surface. Since there are no free currents, the 
magnetic field in the fluid, can be written in terms of a gradient field = h^ — Vcj), 
where h is the magnitude of the applied uniform magnetic field strength in the fluid in 
the vertical direction z and the non-uniform part of the field is described by — V0. The 
magnetisation is given by the linear equation M{Hf) = (/i^ — 1)Hf, where = |Hf| 
and fir =const is the magnetic permeability of the magnetic fluid relative to that of free- 
space, /iQ. Thus the magnetic induction follows from the linear relation B = 

Of course a linear magnetisation law can only be a rough approximation of the real 
M{Hy) measured in the experiment (§3). Despite that, previous analytical studies which 
are based on a linear can successfully predict all regular planforms on the fer- 

rofiuids and their underlying bifurcations (see Gailitis 1977; Friedrichs & Engel 2001; 
Friedrichs 2002; Bohlius et al. 2011). Likewise in our numerical study we aim to catch 
the decisive qualitative features seen in the experiment, and to elucidate whether they 
are connected to homoclinic snaking. 

The magnetic field in the air above the magnetic fluid is given by — Vx and 

the steady state free surface is given by z = C(x, y) with C = 0 corresponding to the fiat 
state. The equations are nondimensionalized using the critical wavenumber kc = \/pg/cF 
for the length scale and y/pgcr for the unit of pressure, where p is the fluid density and a 
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is the surface tension. Furthermore, the amplitude h of the uniform field strength inside 
the magnetic fluid is rescaled using 




hOhr 

2{pr + 1) 


{hr 


l)h, 


(5.1) 


such that the instability occurs at ^ = 1. We note, when comparing equation (5.1) with 
the related (5.5) of the work by Silber & Knobloch (1988), we have rescaled with respect 
to /io- The bifurcation parameter is then 




We note that e depends on fir. For a static magnetic fluid in a vertical field, the dimen¬ 
sionless equations are 

v{p* + z-t(^r- = 0, ^ G {-D, C), (5.2a) 

V^<p = 0, z€{-D,C), (5.26) 

= zeiCD), (5.2c) 


where D is the depth of the fluid and p* is an effective pressure in the magnetic fluid. The 
first of these equations (5.2a) is the stationary Navier-Stokes equation for the magnetic 
fluid, and the other equations (5.26,5.2c) are derived from Maxwell’s equations assuming a 
linear magnetisation law; see Cowley & Rosensweig (1967). We consider a domain above 
the fluid to a height D for computational ease. At the interface z = C, the boundary 
conditions are given by 

P* + ^(Mr - 1)^ (Hf-h)^ = 2k, (5.3a) 

[B-n]+=0, [Hxn]+=0, (5.36) 


where [ • is the value of the jump across the interface of the quantity in the brackets, 
h is the unit normal to the interface, and n is the mean curvature of the surface: 


z-VC _Vh 

(1 + |VC|2)1/2’ 

see Silber & Knobloch (1988). No-flux boundary conditions are imposed on the top and 
bottom of the domain i.e., 0;s(x, —D) = Xz{^^ +T)) = 0 where x = (x,p). Solving (5.2a) 
for the pressure, within a constant, C, at the free-surface yields 

P* = -C+^(Air-l)|HF|2+C', at^ = C. (5.4) 

The pressure can then be eliminated and the boundary conditions written in terms of 
the fields C(a:, p), (/)(a:, p, z), and x(x,p, 2 )), as well as two parameters, fir and e. 

The system of equations possesses a free-energy given by 


[l f 




/ dz + fir 


[2 V 


l-D J 


2 


+CC + hrh{Xz=D — <fz=-D) — a/I + I'^CP “ 1 


dx, 


(5.5) 


made up of magnetic, hydrostatic, surface energies, respectively, with conservation of 
mass, boundary conditions at the top and bottom of the domain, and with the boundary 
condition at the interface 0(x, Q — h({l — fir) = x(x, C) as constraints, where C is 
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the space domain of C(x) and x = {x,y). We note that Gailitis (1977) and Bohlius et al 
(2006) use other energy formulations without the boundary conditions and conservation 
of mass. The constant C is a Lagrange multiplier for the mass constraint: (^dx equal to 

a constant. Taking variations of L with respect to the free-surface ( yields the interface 
equation (5.3a), and taking variations with respect to 0 and y, using the constraint, 
yields the Laplace equations (5.26,5.2c) and the no-flux boundary conditions —D) = 
X^(x, -\-D) = 0. The advantage of this free-energy is that it represents the whole system 
and we can now numerically compute energies taking into account the magnetic fields. 

In order to fix the pressure constant, (7, we assume that the base state (y, 0, = 

(0,0,0) solves the equations corresponding to the mass constraint f^(dx = 0. We choose 
to study the large depth of magnetic fluid case and we set the depth, D, to be 10. For 
l^r = 2, this is yields a depth of 16.42 mm i.e., about 3 times the depth of the fluid in the 
experiment, which is sufficiently large that the domain covering hexagons do not change 
as the depth is increased. We have ignored the effect of the shallow layer in our numerics 
since one would have to solve for an additional magnetic field below the magnetic fluid 
and apply suitable compatibility conditions. 

Carrying out a coordinate transformation to flatten out the free-surface (Twombly & 
Thomas 1983) given by 



and 




y 

D+C 


for C < z < D, 


for -D < z < C 


(5.6) 


(5.7) 


we can perform a Legendre transformation in order to define a Hamiltonian. Arbitrarily 
taking the x-direction, we define the momenta coordinates a = d£/d(px^P = dEjdxx 
and T] = dEjdd^x- Then the Hamiltonian for the system is given by 


D nD 

a(l)xdz + / Pxxdz + r]Cx - ^, (5.8) 

Jo 

on the constraint manifold A4 = {0(x, 0) — h({l — fir) — x(x, 0) = 0, 0a; (x, 0 ) — h(^x{^ ~ 
hr) ~ Xcc(x, 0) = 0}; see Groves et al. (2015) and Groves & Toland (1997) for the similar 
case of water waves. 

The Hamiltonian (5.8) is a constant along solutions in the x-direction and provides a 
wavelength selection principle (see Lloyd et al. 2008, for the same principle applied to 
the Swift-Hohenberg equation). In particular, if we evaluate 1-L along any localized planar 
(hexagon) front connecting to the trivial state C = 0 ftie domain-covering hexagon 
pattern, we find that H = H(0) along the limiting domain-covering hexagon pattern so 
that it belongs to the same level set. Generically, H will only equal H(0) at finitely many 
spatially periodic orbits in the family of spatially periodic patterns and will therefore 
serve as a selection principle that involves only the spatially periodic patterns; see Beck 
et al. (2009) for more details. 


Jo 


6. Numerical results 

The 3D equations (5.26,5.2c) and (5.3a) are solved by first flatting out the free-surface 
using (5.6)-(5.7) so that the computational domains are fixed cuboids. We apply Neu- 
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mann boundary conditions on the four sizes of Q in order to remove the translational 
invariance in x and y that would cause problems with Newton Solvers. This choice of 
boundary conditions on O limits us to cellular hexagons, symmetric planar fronts and 
D 2 symmetric patches involving cellular hexagons computed on the positive quadrant; 
see Lloyd et al. (2008) for more details. To remove the continuous symmetry due to the 
magnetic potentials being defined only up to a constant, we add an extra parameter A 
to the Laplace equations (5.26,5.2c) to be solved along with 0, y and ( and we append 
an extra equation (see Twombly & Thomas 1983) 

x{0,D)+fi4>{0,-D) = 0. (6.1) 

We now have a well-posed boundary value problem with all the continuous symmetries re¬ 
moved allowing us to use Newton methods. We discretise in matlab using fourth-order 
central finite-differences (for planar fronts) or a pseudo-spectral Fourier method (for 
hexagon spikes and patches) in x, ^-directions and Chebyshev pseudo-spectral method 
in ^ with a Newton-GMRES method (Kelley 2003) and secant pseudo-arclength con¬ 
tinuation (Krauskopf et al. 2007). For the Newton-GMRES scheme, one needs a pre¬ 
conditioner. Eor this we pre-compute a sparse incomplete LU-decomposition (with the 
Grout algorithm (Saad 1996)) of the linear system with a flat interface at onset and this 
is fixed for the continuation. The incomplete LU-decomposition is computed with a drop 
tolerance of 1x10“^ and we set a maximum of 20 iterations for GMRES with 40 Arnoldi 
restarts and a tolerance of 10“^. Typical spatial step sizes in x- and ^-directions are 
0.25-0.5 for the Eourier-psuedo spectral method and 20 Chebyshev collocation points in 
2 ) and a fixed arc length step size of 0.01. The numerics are validated by comparing with 
normal form results by Silber & Knobloch (1988). In particular we check that the change 
of sub/supercriticality of rolls and squares occurs at the theoretically predicted values 
of /i^; see the appendix for more details. Eold loci are found by carrying out parameter 
slices and then doing a linear interpolation about the fold values. The Maxwell point is 
computed by appending the Hamiltonian and Energy conditions to the system. 

In the following we present the results obtained with the procedures sketched above. 
We start with the hexagon Maxwell curve (§6.1), which is succeeded by hexagon planar 
fronts (§6.2), hexagon patches (§6.3) and rhomboid patches (§6.4) 

6.1. Hexagon Maxwell curve 

Intuitively, one expects localized patches to exist when the energy of a single hexagon 
cell of the infinitely extended pattern has the same energy as the trivial flat state. This 
occurs at the Maxwell point: a hexagon with a wavelength selected by H has the same 
energy f (hex) = f (0), where hex is a single hexagon cell. Eor close to 1, the hexagon 
Maxwell point was previously computed analytically byEriedrichs & Engel (2001) with 
an energy minimization principle for the wavelength selection. In order to find a Maxwell 
curve relevant for localized hexagon patterns one needs to use the Hamiltonian selec¬ 
tion principle. We find that the Hamiltonian wavelength selection criterion only slightly 
changes the wavelength from that at onset. In figure 7, we trace out the hexagon Maxwell 
curve varying both e and /i^. We find that the Maxwell curve closely follows the fold of 
the domain covering hexagons, in particular for = 2 the Maxwell point is found to be 
e ~ —0.013. Eor close to 1, the height of the hexagons becomes very small making it 
difficult for our numerics to capture since we only compute solutions to a relative error 
of 10“^ and round-off error starts to dominate. We do believe though that the Maxwell 
curve should continue down to = 1. 

At the hexagon Maxwell curve, we expect to see stationary fronts between the flat state 
and domain covering hexagons since we are connecting states with the same energy. 
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Figure 7. Two parameter (e,/ir)-bifurcation diagram showing the existence of the primary 
fold of the domain covering hexagons and their corresponding Maxwell point. 



Figure 8. L2 norm of the displacement of the interface of localized planar hexagon fronts 
undergo homoclinic snaking about the Maxwell point (□) for fir — 2. The (10)-front 
(blue) is snaking in a far larger region of parameter space than the (ll)-front (red). Here 
IKII2 — ]h\ Jo normalized L2 norm of the interface. 







































Homoclinic snaking near the surface instability of a polarizable fluid 13 

Infinitely-many different types of planar fronts are possible depending on the front’s 
interface with respect to the underlying hexagonal lattice. These different fronts can 
be classified by their Bravais-Miller index which we denote by (nin 2 ) where Uj are the 
reciprocals of the intercepts of the line with the lines generated by the lattice vectors 
/j, assigning the reciprocal Uj = 0 if the line does not intersect Rli (see Ashcroft & 
Mermin 1976; Lloyd et al. 2008). It is these different types of fronts that form the sides 
of large hexagon patches. 


6.2. Hexagon planar fronts 

The works of Lloyd et al. (2008); Escaff & Descalzi (2009); Kozyreff & Chapman (2013); 
Dean et al. (2014), suggest that one needs only to consider two main hexagon planar 
fronts, namely the (10) and (11) fronts, in order to understand snaking of fully localised 
patches. We note that in figure 5(f)-(i), we see that the patches are mostly made up 
of (10) fronts, with figure 5(f) also showing the left- and right-hand sides that are (11) 
fronts. Hence, we will concentrate on first understanding the (10) and (11) hexagon planar 
fronts. 

Generically, we expect these hexagon fronts to remain pinned in an open region of 
parameter space where they undergo homoclinic snaking about the Maxwell point (see 
Pomeau 1986; Beck et al. 2009). This is confirmed by our numerical results presented in 
figure 8: For jUr = 2 indeed the (10)- and (ll)-fronts snake about the Maxwell point at 
e ^ —0.013. We find, similar to that reported in other systems (Lloyd et al. 2008; Escaff 
& Descalzi 2009; Kozyreff & Chapman 2013; Lloyd & O’Farrell 2013), that the (10)- front 
snakes over a larger region of parameter space than the (ll)-front. This gives a partial 
explanation as to why mostly patches with (10) sides are observed in our experiments as 
these live in a larger region of parameter space. 

If the computational domain is commensurate with the domain-covering hexagons 
then we expect the snaking branches to terminate near the fold of the domain-covering 
hexagons; see Dawes (2008). In figure 8, we have stopped the snaking branch early and 
indicate that the snaking should continue indefinitely for fronts on the plane. 

6.3. Hexagon patches 

We turn now from the planar fronts to the hexagonal patches as seen in the experiment. 
The numerical results for = 2 displayed in figure 9 show that radial spots bifurcate 
subcritically at onset from the fiat state and undergo a subsequent fold. Along the upper 
branch of the spot there is a De-bifurcation from which a hexagon patch emerges that 
undergoes a series of folds yielding larger and larger patches. We terminate the branch 
early as the snaking branch becomes highly intricate further up the branch (cf. Lloyd 
et al. (2008, Figure 25)). However, we expect this branch to eventually terminate near 
the fold of the domain-covering hexagons. This patch also snakes about around the 
(10)- and (ll)-hexagon fronts (indicated by the shaded regions). We observe that the 
snaking occurs around the hexagon fronts since large patches are effectively made up of 
combinations of hexagon fronts’s sides. The hexagon patch displayed in the central inset 
of figure 9 occurs after the third fold of the hexagon-patch-branch, and comprises seven 
fully developed spikes together with an outer shell of 12 lower humps. It resembles the 
experimental patch presented in figure 1, which is also part of the T 2 = 7.5 s series in 
figure 6. 


6.4. Rhomboid patches 

As seen in the experimental patterns shown in figure 5, the localisation of the patches 
is not necessarily centred around a spike. The patches may equally well organize around a 



Figure 9. Emergence of a fully localized hexagon patch (dashed purple) bifurcating from a 
single spike (green) with = 2. We also plot the snaking regions of the (10)- (blue) and 
(ll)-planar fronts (red). 


valley in-between two spikes. These structures have been called localised rhomboids (Lloyd 
et al. 2008). As shown in figure 10, we are able to find these structures in the same pa¬ 
rameter region as the localised hexagon patches. There are some crucial differences for 
these structures to the localised hexagon patches with a spike at the centre of the patch. 
The first difference is that rhomboid patches do not emerge from a bifurcation of a single 
spike but from a double spike configuration that bifurcates from the fiat state. Figure 10 
explains how patches are grown by first trying to complete a (10)-front side by first adding 
cells in the middle of each of these sides. As the patches get larger, we see the snaking 
becomes more confined to the (10)- and (ll)-front snaking regions. In particular, when 
a patch has only (lO)-front sides, then the snaking diagram makes its largest excursion 
by travelling from the left-most (ll)-front limit to the right-most (10). All other types of 
patches are confined to the (ll)-front region. We note that the existence regions of the 
various configurations are not the same. 

We shall briefly comment on stability of the branch shown in figure 10. It is not possible 
to compute the stability of the branches directly from the numerics since we have not 
written down an equation for the evolution of the fluid velocity. However, we expect the 
branch emanating from e = 0 to initially be unstable since the bifurcation is subcritical. 
We anticipate that the branch then re-stabilises at the first fold to yield a stable two-spike 
patch and between each subsequent fold, the branch will destabilise (due to an eigenvalue 
crossing the imaginary axis at the fold) and then re-stabilise; see (Beck et al 2009). Near 
each fold, we anticipate that there are symmetry-breaking bifurcations leading to non-D 2 
symmetric patches. 

We eventually investigate the effect of changing the depth D on the localised rhomboid 































Figure 10. Emergence of fully localized rhomboid patches bifurcating from the flat state with 
fir = 2. We also plot the snaking regions of the (10)- (blue) and (ll)-planar fronts (red). 



Figure 11. The height of the rhomboid patches for D = 7.5,10, and 15 for = 2 up to the 
second fold. We observe that as we decrease D the height of the patches also decreases but that 
for D ^10 we observe little change in the height. 


patches. In figure 11, we plot the maximum height of the patches for D = 7.5,10, and 15. 
For clarity we are stopping the bifurcation curves before the second fold. As the depth 
increases the snaking diagrams are shifted to the left and the maximum height of the 
patches increases as well. It becomes clear that the existence regions of the localized 
patches are affected by the depth of the fluid. This is believed to be caused by the 
Neumann boundary conditions at the bottom (which insist that the magnetic field decays 
to the applied magnetic field here). In a shallow layer, that constraint also reduces the 
magnetic field within the tip of the spikes, and thus their maximum height. 
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Figure 12. (a) Numerical bifurcation diagram for fir = 3.2 for a two-spike configuration, 
depicted in (b), that has bifurcated off the flat state. The state in (b) was computed for 
B = 11.3 mT. We terminated the branch early. 


7. Comparison between experimental and numerical results 

We now turn to comparing our experimental and numerical results. As stated in §5 our 
model is based on a linear magnetisation law, whereas we presented in §3 the nonlinear 
M{Hf) of the real ferrofluid, with an effective permeability ranging from /ieff = 4.57 at 
B = OmT down to 3.2 at Be = 10.8 mT. Therefore a direct matching of all experimental 
and numerical results can not be expected. For /i^ ^ 1, the linear magnetisation law 
becomes a more accurate approximation of the nonlinear magnetisation law and the 
model predicts homoclinic snaking should persist, albeit with an exponentially small 
width the snaking exists in. 

In order to provide a comparison with the experimental results in §4, we instead chose 
a magnetic permeability = const for our model that captures the qualitative features 
of the experiment. In a first attempt we chose to be the effective permeability at the 
critical point /ieff(77c) = 3.2 (cf. figure 3b) for the model: = 3.2. In figure 12 (a), we 

path-follow the two-spike patch, as depicted in figure 12(b), bifurcating from the fiat 
state to the upper branch. As seen from figure 12 (a), for /i^ = 3.2 localised structures 
exist from B = 10.3 mT to 11.5 mT which is a significantly larger range (about 200%) 
than that observed in figure 6. Moreover the height of the spikes is about 8% larger than 
that observed experimentally. 

In a second attempt, we select a smaller value of magnetic permeability, namely /i^ = 2. 
Here we find that while the critical value of the induction is far from that experimentally 
observed, the width of the bistability region between the fiat state and the domain¬ 
covering hexagons is very close (approximately 0.25 mT in both experiments and numer¬ 
ics); see figure 13. Furthermore, we find from the numerics at the fold of the domain¬ 
covering hexagons the predicted height of the hexagons is 1.2 mm (from figure 4 we find 
the measured height of 1.3 mm) and at onset the height of the stable hexagons is pre¬ 
dicted to be 2.4 mm (and from figure 4 we find the measured height of 2.4 mm) i.e., in 
the bistability region the numerics are within 5% error of that experimentally observed. 
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Figure 13. Maximum height of the domain covering pattern (black line) computed numerically 
for //r = 2 at the critical wavenumber kc ~ 608.94m“^ versus the magnetic induction. The 
maximum height of the rhomboid patches is marked by a brown line. 

Hence, it is possible that phenomenologically = 2 is the appropriate value to take in 
the linear magnetisation law for the numerics for localised hexagon patches. 

The heights of the spikes in the patches are found numerically to be highest at the 
core of the patch and decrease as the spikes get closer to the edge of the patch. As shown 
in Figure 13(a), the maximum height of the rhomboid patches is slightly larger than the 
domain-covering hexagons similar to that seen in the experiments. 

The numerical observation of various co-existing hexagon patches in figure 9, in the 
bistability region of the flat state and the domain-covering hexagons, agrees well with 
the experimental results shown in figure 6 where for the same Hs-values several different 
patches were observed. In particular, we are able to find numerically patches (figure 10 
(a), (b), (c) & (d)) that look to be identical to those seen in figure 5(b), (d), (f) and 
(i). Our numerical methods imposed that the patches should have B 2 symmetry hence 
we were unable to locate any of the other patches seen in figure 5. However, due to the 
existence of a Maxwell point near-by, we believe hexagon spikes/cells can be added to the 
patches to break the D 2 symmetry provided you are in the (ll)-hexagon front snaking 
region. Hence, we have a possible explanation for the existence of all the various types 
of patches observed experimentally. 

Perhaps the largest difference between the numerical and experimental results is that 
the sequence of plateaus for fixed T 2 suggests that one can climb up the snakes for 
increasing H 3 , i.e. the snakes are not vertical as in figure 9 or 10 but slanted, as found 
by Firth et al. (2007) and Dawes (2008) for a conserved quantity. Slanted snaking may 
emerge in the presence of a nonlinear magnetisation law and finite-domain effects and 
could provide an explanation for the emergence of larger and larger patches seen in 
figure 10 . 


8. Conclusion and Outlook 

Utilizing a magnetic pulse sequence we prepare a ferrofiuidic layer to be situated next 
to the unstable branch of the Rosensweig instability. There localized hexagon patches 
emerge spontaneously. The measured sequence of patches indicates homoclinic snaking. 
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This is corroborated by numerics where we investigate the Young-Laplace equation cou¬ 
pled to the Maxwell equations. Homoclinic snaking is demonstrated for two types of 
planar hexagon fronts as well as hexagon and rhomboid patches. The experimentally 
observed sequence of patterns can be interpreted as an imperfect mixture of those pro¬ 
totype scenarios. The numerically unveiled Maxwell point of the hexagons is providing 
an energy argument for those various regular and irregular patches observed since near 
the Maxwell point the energy of a single hexagon spike is the same as the flat state 
and there is no energy penalty in adding/removing hexagon spikes. The width of the 
experimentally observed plateaus is about 0.5% of the control parameter, which is com¬ 
parable to the snaking interval of the numerical prototypes. Moreover, the height of the 
patches in experiment and numerics agrees within 10%. Therefore, we believe that our 
model captures all the most important phenomena of the system. In particular at the 
relative magnetic permeability tends to unity the validity of our model increases. It thus 
unveils an organising centre for the emergence of localised hexagon patches seen in the 
experiment. 

It still remains to explore how various stable configurations of patches connect via 
unstable branches both experimentally (by quasi-statically varying the magnetic induc¬ 
tion B once a patch has been generated) and theoretically. The discrepancy between 
the experimental indications for slanted snaking, and the vertical snaking obtained nu¬ 
merically may be solved by incorporating a nonlinear magnetisation law, finite-size and 
shaiiow iayer effects or by carrying out more experiments to determine the precise extent 
of the piateaus. 

Theoreticaiiy, there are now some very interesting directions one couid take. With 
the existence of an energy functionai for the fuii system that incorporates the boundary 
conditions one couid iook at appiying variationai and spatiai-dynamicai techniques that 
have been deveioped in the water-wave literature; see for instance (Groves & Sun 2008; 
Buffoni et al 2013). 
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Appendix: Validation of the numerical scheme 

We present an overview of the validation results of the numerical scheme. We first piot 
a convergence test for the domain covering hexagons for {e, fir) = (0,2) on the domain 
{x,y) G [—27r, 27r] x [—27r/v^, 27r/v^] and D = 10. In figure 14, we piot the the reiative 
error (reiative to a hexagon spike computed with = Ny = 40) as = Ny = N 
is varied and then as Nz is varied. Since we are using pseudo-spectrai methods, we see 
rapid (geometric) convergence with few mesh points. For our numericai resuits we have 
a reiative error of 10“^ and the foid points are computed to a reiative error of 10“^. 

In order to test that we have correctly captured the nonlinearity in the system, we 
compare our results with the normal form analysis of Silber & Knobloch (1988) shown 
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Figure 14. (a) Bifurcation diagram of the domain covering hexagons for /Xr = 2 on the 
domain {x,y) G [—27r,27r] x [—27r/\/3,27r/\/3] with Nx = Ny = 20,11 = 10 and varying 
Nz. (b) Bifurcation diagram of the domain covering hexagons for = 2 on the domain 
{x,y) G [—27r, 27r] x [—27r/v^, 27r/v^] with Nz = 20 and varying Nx {Ny = Nx). (c) Rela¬ 
tive error on a semi-log scale of the L2 norm of the interface ||C||| for the domain covering 
hexagons on the upper branch at (e,/Xr) = (0, 2) with Nz = 20, H = 10 and varying Nx = Ny 
with respect to the hexagon computed for Nz = 22, Nx = Ny = 22. (d) Relative error on a 
semi-log scale of the L2 norm of the interface HCHi foi' fhe domain covering hexagons on the 
upper branch at (e, /Xr) = (0, 2) with Nx = Ny = 20, H = 10 and varying Nz with respect to the 
hexagon computed for Nz = 22, Nx = Ny = 22. From both panels (c) & (d) we observe rapid 
(geometric) convergence due to the spectral accuracy of the numerical pseudo-spectral methods. 
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Figure 15. Normal form results of Silber & Knobloch (1988) showing the changes in sub- 
and super-criticality of the bifurcations for rolls, squares and up-hexagons i.e hexagons whose 
maximum amplitude is positive. 
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Figure 16. Numerical computation of existence regions of domain-covering rolls, squares and 
hexagons for e < 0. We observe that all three of the domain-covering patterns bifurcate sub- 
critically in the regions where the normal form theory of Silber & Knobloch (1988) predicts. 


in figure 15. Since the analysis is for small amplitude we only compare the bifurcation 
transitions for each of the rolls and squares from sub- to supercritical. Here we find 
excellent agreement with our numerics giving us confidence in our results; shown in 
figure 16 where we plot the existence regions of the rolls, squares and hexagons for e < 0. 
The final validation we did was to find that the Maxwell point computed for = 1.9 
lies in the middle of the snaking region of both the (10)- and (ll)-fronts in figure 8 as 
predicted by (Woods & Champneys 1999; Coullet et al. 2000). This demonstrates that 
the discretisation of (5.26,5.2c) and (5.3a) is consistent with the energy/hamiltonian 
formalisation (5.5) and (5.8) and we are solving the correct equations. 
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